Reduction of the RPA eigenvalue problem and a generalized Cholesky decomposition 

for real-symmetric matrices 
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The particular symmetry of the random-phase-approximation (RPA) matrix has been utilized 
in the past to reduce the RPA eigenvalue problem into a symmetric-matrix problem of half the 
dimension. The condition of positive definiteness of at least one of the matrices A ± B has been 
imposed (where A and B are the submatrices of the RPA matrix) so that, e.g., its square root 
can be found by Cholesky decomposition. In this work, alternative methods are pointed out to 
reduce the RPA problem to a real (not symmetric, in general) problem of half the dimension, with 
the condition of positive definiteness relaxed. One of the methods relies on a generalized Cholesky 
decomposition, valid for non-singular real symmetric matrices. The algorithm is described and a 
corresponding routine in C is given. 

PACS numbers: 21.60. Jz,02.60.Dc,02.10.Yn 



I. INTRODUCTION 



The eigenvalue problem 
approximation (RPA) type, 



of random-phase- 



A B 

-B -A 



= £A 



(1) 



where A and B are real symmetric n x n matrices, is 
frequently encountered in quantum many-body physics. 
In nuclear physics, in particular, the RPA is the most 
widely used method to examine collective excitations of 
nuclei. The RPA is also used to compute the correla- 
tion energy of a system described at zeroth order by the 
Hartree-Fock approximation. Standard first-order RPA, 
relativistic RPA, renormalized RPA, as well as quasi- 
particle and second-order RPA, when formulated in con- 
figuration space, result in equations of the form (JXJ) . The 
properties of the solutions, followingfrom the symmetry 
of the matrix R, are well known [H,H, & 0] • 

The (2n) x (2n) RPA matrix R is not symmetric - 
one reason why for large sizes n it becomes prohibiting 
to solve the RPA problem as-is. It has been possible, 
however, to exploit the particular structure of R, deter- 
mined by the symmetric matrices A and B, in order to 
reduce the RPA problem not only to a symmetric eigen- 
value problem, but at the same time a problem of half the 
dimension, n x n [1, @ • (The generalized RPA problem 
where A is Hermitian and B symmetric can also be re- 
duced Q, but that problem will not be dealt with here.) 
As long as one wishes to avoid complex matrices, that can 
be achieved under the condition that at least one of the 
matrices A + B, A — B be positive definite. The Cholesky 
decomposition of the positive-definite matrix A± B lies 
at the heart of the method presented in Ref. while 
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in Ref. [6( an orthogonal transformation is utilized under 
the same condition. 

It can happen, however, that the positive-definitcncss 
condition of A±B is not met, meaning that the RPA ma- 
trix has imaginary eigenvalues or the stability matrix is 
not positive-definite. In modern relativistic RPA models 
of nuclear response the A± B matrices are known to be 
indefinite, due to the inclusion of states from the Dirac 
sea [8(. Thus the full (2n) x (2n) non-symmetric problem 
is currently solved when using relativistic RPA in config- 
uration space. Other situations from nuclear physics in- 
clude the trivial case of a dipole spurious state appearing 
at imaginary energy, as well as systems which are unsta- 
ble against certain "excitations" . An exotic nucleus with 
different enough proton- and neutron- Fermi energies can 
be found unstable against configurations of isospin T = 1 
and certain angular momentum and parity J* . Also in 
the ncighborgood of phase transitions the stability matrix 
can have negative eigenvalues when the residual interac- 
tion is attractive and large enough. In such cases the 
(in)definiteness of the matrices may not be known before 
solving the RPA problem. 

In this work alternative procedures are proposed for 
reducing the size of the RPA problem, still involving real 
matrices, where the condition of positive-definitness is 
relaxed. The result is a real, non-symmetric (in general) 
problem of half the original dimension. One procedure re- 
quires no matrix decomposition and is recommended for 
problems where the A± B matrix is expected from the 
outset to not be positive-definite. Another method offers 
the possibility to detect and solve a symmetric problem 
if the matrix turns out to be positive-definite. The latter 
involves a more general Cholesky-like decomposition of 
a real symmetric matrix. The way to perform the de- 
composition in practice is outlined and a C routine is 
also provided. Compared with the usual Cholesky de- 
composition, the additional computational and storage 
effort is minimal. Other reduction methods can be de- 
vised based on different matrix decompositions, as will 
be demonstrated. 
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Modern computers perform quite well in solving the 
standard first-order RPA in nuclei, without reduction. 
Still, the need to reduce the RPA problem and thus ac- 
celerate its solution becomes rather pressing when, e.g., 
one solves the RPA equations iteratively (as in renormal- 
ized RPA) or one wants to evaluate the RPA correlation 
energy of a nucleus. In such cases many RPA problems 
corresponding to different JT quantum numbers must 
be solved (perhaps more than once) before a solution is 
reached. Hence a fast technique is of great importance. 

Reduction methods are of interest also when one needs 
to evaluate only the lowest positive eigenvalues of very 
large RPA matrices, using, e.g., the Lanczos technique - 
see Ref. [§] for a related strategy. 

Next it is briefly reviewed how the RPA problem can 
be reduced using the Cholesky decomposition. Then al- 
ternative methods and the modified Cholesky decompo- 
sition and its algorithm are presented, a corresponding 
routine in C is given and comments follow before con- 
cluding. 



Once the eigenvectors R\ have been evaluated, the vec- 
tors Xx and Y\ can be recovered. For real and positive 
eigenvalues (and real eigenvectors), 



i 



T\-l 



L]R\, 



L]R X . 



(8) 
(9) 



It is easily verified that the desired normalization condi- 
tion as well as the orthogonality condition between dif- 
ferent eigenvectors 



XiX ll -Y^Y li = S 



(10) 



are satisfied. 

If the matrix A+B is not positive definite, but A—B is, 
one can write A — B = LL T and proceed in an analogous 
way by rewriting the first equation of the system ([3]) in 
the form §5$ with 



II. REDUCTION OF THE RPA PROBLEM 
USING CHOLESKY DECOMPOSITION 

The eigenvalue problem of eq. {1} leads in a straight- 
forward manner to the system of equations 



H = L T (A + B)L 



-1/2 T T 



L 1 (X, 



Y x 



From R\ the vectors X\ and Y\ can be recovered, 
1, 1 



X x = -\ 



L + t/£~x(L t )~ 1 ]Rx, 



(11) 



(12) 



(A-B)(Xx-Yx) =ex(Xx + Y x ) 
(A + B){Xx+Yx)=e x (Xx-Yx) 



(A + B)(A-B)(Xx-Y x 
(A-B)(A + B)(Xx+Yx 



= e\{Xx-Yx) 
= el(Xx + Yx) 



(2) 



(3) 



The matrices A ± B are real and symmetric. Let us 
assume that A + B is positive-definite. Then it can be 
factorized as 



-x^r^Rx. 



(13) 



For more details see Ref. @. 

Among the virtues of this method of solving the RPA 
problem are that the Cholesky decomposition is fast and 
efficient numerically [l(| and that one only has to deal 
with triangular and symmetric matrices, which simplifies 
the numerical realization of the solution. 



A + B = LL , 



(4) 



where L a lower-triangular real matrix and L T its trans- 
pose. This is the square-root or Cholesky decomposition 
of the matrix. Then the second equation of the system 
([3]), premultiplied with L T , can be written as 



where 



HRx 



H = L T (A — B)L 



a real symmetric matrix and 



Rx 



^ 2 L T (Xx + Yx) 



(5) 



(6) 



(7) 



its orthonormalized eigenvectors with eigenvalues e\. 

The original RPA problem has been reduced to a sym- 
metric problem of half the dimension. For e\ > the 
eigenvalues of the original problem are real, ex = ± 



III. GENERALIZED METHODS 

Unfortunately, if neither A + B nor A — B is posi- 
tive definite, the above method fails and one has to ei- 
ther use general routines to solve the full non-symmetric 
(2n) x (2ri) RPA problem or use complex n x n matrices 
(see Ref. [6[ ; alternatively, one can define a Cholesky de- 
composition LL T where imaginary values are allowed in 

One can still reduce the size of the problem using other 
decompositions, valid for generic real symmetric matrices 
[DI EH) EH or > m fact, no decomposition at all, as will be 
shown. For the moment, the discussion is kept general. 
Let us consider the second equation of the system 
One can of course start with the first one and proceed 
in a completely analogous way. Suppose we can factorize 
the matrix (A + B) as 



A + B = CDE, 



(14) 
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where C, D, E are n x n matrices and E has an inverse. 
(Any one of them may be the identity matrix I). The 
equation in question can be written as 



HR\ 



with 



H = E(A - B)CD , R x 



(15) 



E(X X + Y X ). (16) 



Henceforth we consider the solutions with real and posi- 
tive e\. Real eigenvalues allow us to assume real vectors 
R x . The real vectors 

Xx = k^E- 1 + -^CD]R X , (17) 



Y x 



form solutions of the RPA problem, eq. |T]), as can be 
easily verified. They obey the orthonormalization condi- 
tion 



kV^E- 1 - -^—CD]R X (18) 
1 V £ x 



which is equivalent to 

rKe-Ycdr^ = ±s Xlt . 

Note that 



(19) 



(20) 



(E~ 1 ) T CD = D T C T E- 1 = {E^fiA + B)E~ 1 . (21) 



by writing 
, then writ- 



Expressions (fT7|) . (JTHJ) can be obtained 
{X x oxY x } = l[(X x +Y x ){+or-}(X x -Y x 
ing (X x -Y x )=e x 1 (A + B){X X + Y x ) (see eq. ©) and 
finally using R x = ^/e x E^ 1 (X x + Y x ) and the decompo- 
sition of (A + B). Eq. ([20]) is easily obtained if one sub- 
stitutes X and Y in eq. (Till|) with expressions (jTTJ) and 
(JTSJ) . Thus, the existence of the orthonormalizable real 
solutions (X Xtft , Y x>fl ) for the RPA matrix and of the in- 
verse of E guarantee the orthogonality condition eq. (j2"D|) 
for the corresponding eigenvectors of H (which has been 
verified by means of numerical examples for the special 
cases presented in the next two subsections). 

The normalization condition does not hold automati- 
cally. Therefore, either the eigenvectors R x have to be 
renormalized according to eq. (j20|) (for A = /i) after solv- 
ing the eigenvalue problem (fTS"]) . or the X x and Y x vectors 
have to be renormalized according to eq. (|19p , after they 
are calculated from R x and eqs. (jTTJ), (|18p. Regarding 
the r.h.s. of eq. (fT9| . note that the norm of the RPA 
eigenvectors - as defined by the l.h.s. - need not be posi- 
tive, unless the stability matrix is positive definite [2j. If 
it is not, eigenvectors with positive eigenvalues may exist 
which arc normalizable to —1, instead of 1. Such a case 
can be recognized before normalization by evaluating the 
l.h.s. of expression (Ti~9"l) or (f2T)]) and checking the sign of 
the result. 



In Ref. @| an orthogonal transformation was utilized, 
A±B = CDC T , where C orthogonal and D diagonal. (It 
was applied in a different way so as to obtain a symmet- 
ric problem - real or complex.) In the method described 
previously based on the usual Cholesky decomposition, 
one has C = L, D = I and E = L T and the generic 
equations (|16p - (l2"01) simplify into the ones presented ear- 
lier and in Ref. [5j . Indeed, a procedure such as outlined 
above would make less sense to apply if the matrices C, 
D, E did not have special properties which simplify the 
algebra and numerics involved. 



A. No decomposition 

One may chose simply to solve one of the n X n 
eigenvalue problems of eq. ^ without any decomposi- 
tion. For example, one can solve the second one for 



.-1/2 



(X x + Y x ) and apply the equations given 



above for 



D = E 



C = A + B. 



(22) 



Some matrix operations are thus saved, but the re- 
sulting eigenvalue problem is not symmetric. The no- 
decomposition strategy should be preferable when the 
matrices are expected to not be positive-definite. 



B. Generalized Cholesky decomposition 

A decomposition strategy like the one discussed next 
offers the additional possibility to detect and solve a real- 
symmetric problem, if the decomposed matrix turns out 
to be positive-definite. It can be shown that a non- 
singular real symmetric matrix F with non-singular lead- 
ing submatrices, with or without negative eigenvalues, 
can be factorized in the form 



F = LDL , 



(23) 



where L is again a lower-triangular matrix and D = 
diag{c?i} is a diagonal matrix whose diagonal elements 
are equal to 1 or —1. Indeed, the factorization (f2"3")) fol- 
lows from the "LZ3L" decomposition of linear algebra, 
let us write it as 



F = L'D'L' T , 



(24) 



where L' is a unit-triangular matrix and D' = diagjc^} 
a general diagonal matrix [l2|, [H[ . If we write D' as 

D' = [dmg{^\}} [diag{sgnda] [diagly^}] 

and define L and D as 

L = L' [diag{y^]}] , D = diag{sgn<}, 

eq. (|23p is obtained. The existence of the LL LDV de- 
composition (eq. (|24p ) for a non-singular real symmetric 
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matrix F with non-singular leading submatrices is known 
from linear algebra [T2, EH • 

Setting C — L and E = L T in eqs. fH)] ) -(f2Ti ]) we obtain 
the n x n real, non-symmetric (in general) eigenvalue 
problem of the form (|15|) with 

H = L T (A - B)LD , R x = -Ll t (X x + Y x ). (25) 

At this point one may check whether D = I (meaning 
that F is positive-definite) , in which case H is symmetric 
and optimized routines can be used to solve its eigenvalue 
problem. 

Applying eqs. (jTTJ) (fT5|) we find that for real and pos- 
itive eigenvalues the solutions of the RPA equations are 
given by the vectors 

A a - kv^i^r 1 + ~^LD}R X , (26) 

Yx = -U/iX^r 1 - ~^=LD]R X . (27) 

The normalization condition (|20|) is simplified since 
(E-YCD = D. 

The "generalized" Cholesky decomposition defined by 
eq. (|23p can be realized numerically almost as efficiently 
as the usual one. A simple code in C is given later on. 
The underlying algorithm (a revision of a similar one 
which appears in Ref. (HI) resembles closely the usual 
Cholesky algorithm and can be described as follows. 

If we write out eq. (|23|) in components (F = [fij], 
D = [dij] with dij — diSij, L — [hj]), we have 

fij = ^ dkkkljk, (28) 

k<i,j 

from which we readily obtain, setting i = j and i < j 
respectively, 



V k<i 



Iji — j~{fj% — ~^^dkhkljk)- (30) 

Hi 



We start the solution with In — \fd\f\i = \f\fn\] if 
/n < 0, we have d\ = — 1, otherwise d\ = 1. The off- 
diagonal elements lj\ can now be evaluated. After fin- 
ishing with the first column of L, we continue with Z22 ■ 
Similarly, the sign of di will be determined by the sign 
of (/22 — cMii)- All quantities needed to calculate Z j2 are 
now known. In short, wc apply eqs. (|30[) in the order 
% = 1, 2, . . . , n, each time starting with the diagonal el- 
ement la and proceeding with the off-diagonal elements 
in the same column. 



IV. A ROUTINE IN C 

A simple routine example is given below in C language, 
based on the choldc routine of fioj ] . § 2.9 and the al- 
gorithm described previously. Given the real-symmetric 
n x n matrix F (f [1 . . . n] [1 . . .n] ), this routine con- 
structs its decomposition F = LDL T , where D is a di- 
agonal matrix with elements equal to 1 or —1 along the 
diagonal and L a lower-triangular matrix. Only com- 
ponents F^ with j > i need to be referenced, since F 
is symmetric. This allows the elements of L lying below 
the diagonal to be stored in the corresponding elements of 
f . Additional arrays p [1 . . . n] and d [1 . . . n] are used to 
store the diagonals of L and D respectively. One can also 
use the routine to test if the matrix F is positive-definite, 
by looking for —1 elements in D. The return value can, 
e.g., be an integer (instead of void) equal to the number 
of negative d [i] values. For matrices which are expected 
to have mostly positive eigenvalues, in which case most 
d[i] 's will be equal to +1, it should be advantageous to 
only store the i-values for which d[i]= —1, instead of 
the whole array d, especially if the matrix dimension is 
large. 

float f [N] [N] ; // N=n+1 
float p[N] ; 
int d [N] ; 

void choldmod(int n) 
{ 

int i,j,k; 
float sum,aux; 

for (i=l;i<=n;i++) 
{ 

d[i]=l; 

for (j=i; j<=n; j++) 
{ 

sum=0 . ; 

for(k=l; k<i; k++) 

sum += d [k] *f [i] [k] *f [ j ] [k] ; 
aux = f [i] [j]-sum; 
if (i==j) 
{ 

if (aux<=0.0) d[i]=-l; 
p [i] =sqrt (d [i] *aux) ; 

} 

else 

f [j] [i]=d[i]*aux/p[i] ; 

} 

> 

} 

The inner loop, consisting of two multiplies (one of 
them with int) and a sum, is executed (n 3 — n)/6 rt 3 /6 
times. There are also n square roots and about n 2 /2 di- 
vides. On an Intel Pentium 4 machine (gec compilation) 
an arbitrary 1000 x 1000 symmetric matrix was factorized 
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in about 2 seconds and a 10000 x 10000 matrix in about 
40 minutes. The algorithm of the usual Cholesky decom- 
position contains n(n 2 — l)/6 ~ n 3 /6 similar inner loops 
(with one multiply and one subtract each) , n square roots 
and about n 2 /2 divides [Toj ] . A similar amount of effort 
is required for a usual LDL decomposition, whereas the 
orthogonal transformation proceeding via a Householder 
reduction and a QL decomposition requires more effort, 
namely about 2n 3 /3 + 30n 2 operations [lOj . As holds for 
the usual LDL factorization, stability of the present al- 
gorithm is not guaranteed for indefinite matrices. Since, 
however, RPA matrices tend to have their largest ele- 
ments (absolute values) along the diagonal, they should 
be mostly well behaved. 

V. COMMENTS 

In most of the above it is assumed that there are no 
singularities. In actual applications it is quite improbable 
to obtain eigenvalues so close to zero that bad behaviour 
occurs. Even then, one can escape the pitfall by slightly 
scaling the residual interaction entering the RPA equa- 
tion and repeating the calculation. 

The modified eigenvalue problem defined by eqs. (|15[) 
and (|25[) can be transformed into a complex symmetric 
problem. One has to multiply from the left both sides 
of eq. lfT5| with the square-root matrix of D, i.e., with 
the diagonal matrix D with diagonal elements di = \fd~i 
equal to 1 or i, so that D = D 2 . The i— th element of 
the eigenvector of the new, symmetric (but complex, if 
D contains negative elements) matrix DL T (A — B)LD, 
is equal to the i— th element of the corresponding eigen- 



vector of H, times the i— th diagonal element of D. That 
the RPA problem can be reduced in a complex symmetric 
problem was shown in a different way in Ref. [6(. 

The modified Cholesky decomposition defined here is 
not available in packages like LAPACK, appropriate for 
very large matrices, contrary to the usual Cholesky and 
other decompositions. Therefore, for very large matrices 
one may have to chose another reduction procedure. 



VI. CONCLUSION 

In conclusion, it has been demonstrated that there ex- 
ist methods to reduce a real RPA eigenvalue problem to a 
real non-symmetric problem of half the dimension with- 
out demanding that one of the matrices A + B, A — B 
be positive definite. Reduction can be achieved with or 
without matrix decomposition. We have worked out a 
method based on a generalized Cholesky decomposition, 
which is no more involved numerically than the one re- 
lying on the usual square-root decomposition of positive- 
definite real-symmetric matrices. The result is in general 
a real non-symmetric (or a complex-symmetric) eigen- 
value problem of half the dimension. If one of the ma- 
trices turns out to be positive definite, a real-symmetric 
problem is obtained instead, without additional effort. 
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